################################################################################
###
### Filename: 	Map-Racial-Liberalism.R
### Author: 	  David A. Bateman
### Function: 	Reproduces maps of support for abolitionism and for free Black
###				      rights in Virginia, North Carolina, Tennessee, and Mississippi.
###             These maps are included only in the online appendix.
###             Parts of this code are based on slavery-maps-replication.R, 
###             submitted by Acharya, Blackwell, and Sen to the Journal of 
###             Politics dataverse to reproduce figures and tables in Acharya, 
###             Blackwell, and Sen (2016). 
###	Last update: January 1, 2023
###
################################################################################

rm(list = ls())

library(readxl)
library(dismo)
library(ggmap) 
library(raster)
library(rgdal)
library(rgeos)
library(elevatr)
library(plyr)
library(reshape)
library(ineq)
library(maps)
library(maptools)
library(foreign)
library(RColorBrewer)
library(classInt)
library(shapefiles)
library(ggplot2)
data(state.fips)
state.fips <- unique(state.fips[,c("fips","abb")])
state.fips$abb <- as.character(state.fips$abb)
state.fips <- rbind(state.fips, c(2, "AK"))
state.fips <- rbind(state.fips, c(15, "HI"))
rownames(state.fips) <- state.fips$abb
fips.state <- state.fips
rownames(fips.state) <- fips.state$fips
data(county.fips)

dir <-paste("../Deeper Roots/", sep="")
setwd(dir)

## Historical county data available from http://publications.newberry.org/ahcbp/
histcounties <- readOGR("Dataverse/Shapefiles/US_AtlasHCB_Counties_Gen001/US_HistCounties_Gen001_Shapefile", layer = "US_HistCounties_Gen001")
histcounties$fips <- as.numeric(as.character(histcounties$FIPS))
histstates <- readOGR("Dataverse/Shapefiles/US_AtlasHCB_StateTerr_Gen0001/US_HistStateTerr_Gen0001_Shapefile", layer = "US_HistStateTerr_Gen0001")
ss2000 <- which(histstates$START_N <= 20001231 & histstates$END_N >= 20001231)
states2000 <- histstates[ss2000,]
states2000 <- states2000[!(states2000$ABBR_NAME %in% c("AK", "HI")),]

## Restrict to those states included in the analyses of voting on
## rights for free persons of color. 
voting.states <- c("Virginia", "Tennessee", "West Virginia", "North Carolina", "Mississippi")
## Restrict to those states which voted on considering or preventing abolition 
abolition.states <- c("Virginia", "Tennessee", "West Virginia")

## Start with rights for free persons of color.
south <- states2000$ABBR_NAME %in% c("VA", "TN", "WV", "NC", "MS")
states2000$south <- south
the.south <- unionSpatialPolygons(states2000, south)

## These data sources are originally from Jeremy Atack. The links and code are 
## from interpolation.R (Acharya, Blackwell, and Sen (2016)).
## https://my.vanderbilt.edu/jeremyatack/data-downloads/
rivers <- readOGR("Dataverse/Shapefiles/SteamboatNavigatedRiverContigUSAprojection", layer = "SteamboatNavigatedRiverContigUSAprojection")
rivers <- rivers[rivers$START < 1830 & rivers$END > 1840,]
canals <- readOGR("Dataverse/Shapefiles/19thC_Canals_March2017", layer = "19thC_Canals_March2017")
canals <- canals[which(canals$OPENED < 1830 & canals$CLOSED > 1830),]

## 1830 Census
## Data comes from ICPSR 35206
## http://www.icpsr.umich.edu/icpsrweb/ICPSR/studies/35206
census1830 <- read.dta("Dataverse/Census/census1830.dta")
vartable <- cbind(names(census1830), attr(census1830, "var.labels"))
cvars <- c("totpop", "stot", "fips")
census1830 <- census1830[census1830$level == 1, ]
census1830 <- census1830[, cvars]
census1830 <- census1830[!is.na(census1830$fips),]
census1830$fips <- floor(census1830$fips)
census1830 <- cast(melt(census1830, id = "fips"), fips ~ variable, sum)
census1830$sprop <- 100*(census1830$stot/census1830$totpop)

counties.map <- spTransform(histcounties, CRS(proj4string(rivers)))
canals <- spTransform(canals, CRS(proj4string(rivers)))
the.south <- spTransform(the.south, CRS(proj4string(rivers)))

counties.map$fips[which(counties.map$NAME == "CARROLL (ext)")] <- 22035

cc1830 <- which(counties.map$START_N <= 18290101 & counties.map$END_N >= 18290101)
counties1830 <- counties.map[cc1830,]
counties1830@data$order <- 1:nrow(counties1830@data)
counties1830@data <- merge(counties1830@data,census1830, by=c("fips"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties1830@data<-counties1830@data[order(counties1830@data$order),]

### Free persons of color rights counties
votes <- read.csv("Dataverse/Files/County-Votes-Black-Rights.csv")
liberal.counties <- merge(counties1830, votes, by=c("fips"), all.x=TRUE, all.y=FALSE, sort=FALSE)
liberal.counties <- liberal.counties$new>=.5
counties1830$liberal.counties <- liberal.counties
liberal.counties <- unionSpatialPolygons(counties1830, liberal.counties)
liberal.counties <- spTransform(liberal.counties, CRS(proj4string(rivers)))

### Abolition counties
abolition <- read.csv("Dataverse/Files/Virginia-TN-Abolition.csv")
abol.counties <- merge(counties1830, abolition, by=c("fips"), all.x=TRUE, all.y=FALSE, sort=FALSE)
abol.counties <- abol.counties$abolition>=.5
counties1830$abol.counties <- abol.counties
abol.counties <- unionSpatialPolygons(counties1830, abol.counties)
abol.counties <- spTransform(abol.counties, CRS(proj4string(rivers)))

## Create a heat map function
plot.heat1 <- function(counties.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- color.palette1(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map, border=NA, lwd=1,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat), add=FALSE)
  
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
}

color.palette1 =  colorRampPalette(c("white","grey","black"),space="rgb", bias=1)

################################################################################
### Figure A16: Votes on abolition, dashed counties’ representatives votes   ###
###             in favor.                                                    ###
################################################################################
label<-paste("Dataverse/Figures/Figure-A16.pdf", sep="")
pdf(label, width=8, height=6, bg="transparent", family="sans" )
par(mfrow=c(1,1),oma=c(1,0,1,0), mar=c(0,0,2,0))

plot.heat1(counties1830[counties1830$STATE_TERR %in% abolition.states,],z="sprop",breaks=c(0:100), plot.legend=FALSE)

plot(rivers[rivers$END > 1830,], col = "dodgerblue", lwd = 1, add = TRUE)
plot(canals, col = "blue", lwd = 1, add = TRUE)
plot(abol.counties[2,], density=c(10), lty=2, add=TRUE)
dev.off()

################################################################################
### Figure A17: Votes on civil disabilities for free African Americans,      ###
###             dashed counties’ representatives voted against.              ###
################################################################################
label<-paste("Dataverse/Figures/Figure-A17.pdf", sep="")
pdf(label, width=8, height=6, bg="transparent", family="sans" )
par(mfrow=c(1,1),oma=c(1,0,1,0), mar=c(0,0,2,0))

plot.heat1(counties1830[counties1830$STATE_TERR %in% voting.states,],z="sprop",breaks=c(0:100), plot.legend=FALSE)

plot(rivers[rivers$END > 1830,], col = "dodgerblue", lwd = 1, add = TRUE)
plot(canals, col = "blue", lwd = 1, add = TRUE)
plot(liberal.counties[1,], density=c(10), lty=2, add=TRUE)
dev.off()
